This script is a focused analysis of mesenchymal to epithelial signaling in the BEFM project data, re-run on a deeper standardized sampling (5000 cells per sample), looking specifically at BEF12, BEF14, BEF15, BEFM1 and BEFM2 (per Allie’s request 2023-11-08.)

We will isolate and then look at Mesenchymal Influence, the Epithelial Microenvironment, and Mesencymal-Epithelial Signaling. These three indepednet lenses should allow us to clearly understand what is going on.

This script is a working draft of what might evolve to become a useful standardized framework to look at the full dataset, later. This workflow is experimental for now and is in active revision in collaboration with Allie.

Setup Packages

Here we load our required packages for this analysis.

# Set WD
setwd("/Users/msbr/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/.shortcut-targets-by-id/1VLqBlyzO-Qad5O2kbwXkBRh_1cQho39t/Engineered Lung Paper/Global_Connectomics")

# Require packages
require(Seurat)
require(NICHES)
require(ggplot2)
require(cowplot)
require(dplyr)
require(knitr)

Load NICHES data from Part 1

These data have not been filtered, per Part 2 decisions. This seems to be OK for this project, and generally preferable, because some of our measurements, particularly from BCL5, are low-information and we do not want to lose these data points. So, we load from the last save-point, which contains data from a (maximally) deep sampling of the same number of cells from each sample.

load("~/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/.shortcut-targets-by-id/1VLqBlyzO-Qad5O2kbwXkBRh_1cQho39t/Engineered Lung Paper/Global_Connectomics/global.connectomics.2023-11-11.Robj")

Subset out the just specific samples of interest, per Allie’s request:

Here we limit all of our downstream analysis to information that came exclusively from the following 5 samples of interest, which will allow us to nicely isolate the effect of macrophages on mesenchymal to epithelial signaling.

samples.of.interest <- c('BEF12','BEF14','BEF15','BEFM1','BEFM2')
# we will use this set of samples for downstream calculations in this workflow

Define colors

This is a temporary color palette, and we should consder updating for publication.

load("~/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/.shortcut-targets-by-id/1VLqBlyzO-Qad5O2kbwXkBRh_1cQho39t/Engineered Lung Paper/Allie's Object Filtration/Final Objects/all.color_palettes_2.R")

## For now:
color_pals$type_colors <- c('#A40606','#9CFFFA','#B0DB43','#9C528B','#2F6690',
                                '#946846','#F1C40F','green','#0F0326','#E65F5C','#14591D','#726DA8',
                                'yellow','purple','blue','red','orange','darkgrey','magenta')
names(color_pals$type_colors) <- unique(global.connectomics$CellToCell$alra$SendingType)

Isolate Mesenchymal Influence

This is just the Mesenchymal Cell-To-System data, which incorporate information about both cell-signaling character and cell population presence / absence / representation.

temp <- global.connectomics$CellToSystem$alra # let's just look at the imputed data, for now, for ease
Idents(temp) <- temp$class
table(Idents(temp))
## 
##  Epithelium Endothelium  Mesenchyme      Immune 
##       32888        8164       31702        5071
mes.inf <- subset(temp, idents = 'Mesenchyme')
Idents(mes.inf) <- mes.inf$orig.ident
table(Idents(mes.inf))
## 
##  BEF1  BEF2  BEF3 BEF12 BEF14 BEF15 BEFM1 BEFM2 BEFM4 BEFM5 BEFM6  FB13  FB14 
##  2514   742  1520  4497  4037  4066   986  3250    85    82   143  4878  4902
mes.inf <- subset(mes.inf,idents = c(samples.of.interest))
table(Idents(mes.inf))
## 
## BEF12 BEF14 BEF15 BEFM1 BEFM2 
##  4497  4037  4066   986  3250

Isolate Epithelial Microenvironment

This is just the Epithelial System-To-Cell data, which incorporates information about all signaling landing on Epithelium, regardless of source. It changes, like above, both with changing cellular phenotype and with altered cell population representation.

temp <- global.connectomics$SystemToCell$alra
Idents(temp) <- temp$class
table(Idents(temp))
## 
##  Epithelium Endothelium  Mesenchyme      Immune 
##       32888        8164       31702        5071
epi.mic <- subset(temp, idents = 'Epithelium')
Idents(epi.mic) <- epi.mic$orig.ident
table(Idents(epi.mic))
## 
## BC1P3 BC1P6  BCL5 BCEC2  BEF1  BEF2  BEF3 BEF12 BEF14 BEF15 BEFM1 BEFM2 BEFM4 
##  5000  2887  5000  4967   751   451   632   350   806   781  1502   788   752 
## BEFM5 BEFM6 
##  3875  4346
epi.mic <- subset(epi.mic,idents = c(samples.of.interest))
table(Idents(epi.mic))
## 
## BEF12 BEF14 BEF15 BEFM1 BEFM2 
##   350   806   781  1502   788

Isolate Mesenchymal to Epithelial Signaling

This just looks at the Mesenchymal-To-Epithelium Cell-To-Cell communication. It predominantly captures signal due to altered cellular character. It does not appropriately capture altered cell-signaling cues due to altered cell population representation, but it MAY capture signal present in this cross DUE to the addition of unseen cells outside of this cross (i.e., it could change with/without the presence of macrophages, even though there will be no data FROM macrophages in this object.)

temp <- global.connectomics$CellToCell$alra
Idents(temp) <- temp$class.Joint
table(Idents(temp))
## 
##   Epithelium - Epithelium  Epithelium - Endothelium  Endothelium - Epithelium 
##                     87508                     12249                     12249 
## Endothelium - Endothelium   Epithelium - Mesenchyme  Endothelium - Mesenchyme 
##                     26980                     17089                      8196 
##   Mesenchyme - Epithelium  Mesenchyme - Endothelium   Mesenchyme - Mesenchyme 
##                     17089                      8196                     55356 
##       Epithelium - Immune      Endothelium - Immune       Mesenchyme - Immune 
##                      4343                      3040                      2050 
##       Immune - Epithelium      Immune - Endothelium       Immune - Mesenchyme 
##                      4343                      3040                      2050 
##           Immune - Immune 
##                      9763
mes.epi <- subset(temp,idents = 'Mesenchyme - Epithelium')
Idents(mes.epi) <- mes.epi$orig.ident
table(Idents(mes.epi))
## 
##  BEF1  BEF2  BEF3 BEF12 BEF14 BEF15 BEFM1 BEFM2 BEFM4 BEFM5 BEFM6 
##  2212  1071  1376  1386  2578  2540  2511  2134   360   335   586
mes.epi <- subset(mes.epi,idents = c(samples.of.interest))
table(Idents(mes.epi))
## 
## BEF12 BEF14 BEF15 BEFM1 BEFM2 
##  1386  2578  2540  2511  2134

Ok, now we have isolated our data of interest. Next, let’s visualize each data set to get a first sense of ‘where’ there might be signal in the data that we are interested in exploring more deeply:

Mesenchymal Influence

mes.inf <- ScaleData(mes.inf)
mes.inf <- FindVariableFeatures(mes.inf)
mes.inf <- RunPCA(mes.inf,npcs = 100)
ElbowPlot(mes.inf,ndims = 100)

mes.inf <- RunUMAP(mes.inf,dims = 1:25)

p1 <- DimPlot(mes.inf,group.by = 'Dataset2',cols=color_pals$dataset_colors,raster=F,shuffle = T)
p2 <- DimPlot(mes.inf,group.by = 'orig.ident',cols = color_pals$sample_colors,raster=F,shuffle = T,label=T)
p3 <- DimPlot(mes.inf,group.by = 'class',cols = color_pals$class_colors,raster=F,shuffle = T)
p4 <- DimPlot(mes.inf,group.by = 'SendingType',cols = color_pals$type_colors,raster=F,shuffle = T)

plot_grid(p1,p2,p3,p4)

pdf(file='MesenchymeInfluence_Subset.UMAPS.pdf',width=10,height=8)
print(plot_grid(p1,p2,p3,p4))
dev.off()
## quartz_off_screen 
##                 2

I think this look pretty good. Significant tissue-tissue separation, which may be a consequence of real biology or may be batch, we do not know yet. We could try integrating these later if we would like, not sure if it would help us or obscure interesting patterns. Let’s avoid for the moment, but can revisit this decision.

Note that the above plot shows how fibroblast influence is different between samples, but considers all signals, not just those landing on epithelial cells. A portion of the variance shown here is likely due to changes in signaling landing on epitehlial cells. We do not know yet.

Epithelial Microenvironment

epi.mic <- ScaleData(epi.mic)
epi.mic <- FindVariableFeatures(epi.mic)
epi.mic <- RunPCA(epi.mic,npcs = 100)
ElbowPlot(epi.mic,ndims = 100)

epi.mic <- RunUMAP(epi.mic,dims = 1:25)

p1 <- DimPlot(epi.mic,group.by = 'Dataset2',cols=color_pals$dataset_colors,raster=F,shuffle = T)
p2 <- DimPlot(epi.mic,group.by = 'orig.ident',cols = color_pals$sample_colors,raster=F,shuffle = T,label=T)
p3 <- DimPlot(epi.mic,group.by = 'class',cols = color_pals$class_colors,raster=F,shuffle = T)
p4 <- DimPlot(epi.mic,group.by = 'ReceivingType',cols = color_pals$type_colors,raster=F,shuffle = T)

plot_grid(p1,p2,p3,p4)

pdf(file='EpithelialMicroenviornment_Subset.UMAPS.pdf',width=10,height=8)
print(plot_grid(p1,p2,p3,p4))
dev.off()
## quartz_off_screen 
##                 2

This is actually pretty intresting and I hven’t noticed these patterns before. The tissues group more together here. The similarity between Mono + Tri_E vs. Co + Tri_L + Quad is interesting and I suspect relevant.

BEFM1 is clearly different from the rest, suggesting a relatively unique epithelial microenvironment. This aligns with our biological intuition, given the histology and Allie’s analysis.

Mesenchymal-Epithelial Communication

mes.epi <- ScaleData(mes.epi)
mes.epi <- FindVariableFeatures(mes.epi)
mes.epi <- RunPCA(mes.epi,npcs = 100)
ElbowPlot(mes.epi,ndims = 100)

mes.epi <- RunUMAP(mes.epi,dims = 1:20)

p1 <- DimPlot(mes.epi,group.by = 'Dataset2.Sending',cols=color_pals$dataset_colors,raster=F,shuffle = T)
p2 <- DimPlot(mes.epi,group.by = 'orig.ident',cols = color_pals$sample_colors,raster=F,shuffle = T)
p3 <- DimPlot(mes.epi,group.by = 'class.Sending',cols = color_pals$class_colors,raster=F,shuffle = T)
p4 <- DimPlot(mes.epi,group.by = 'class.Receiving',cols = color_pals$class_colors,raster=F,shuffle = T)
p5 <- DimPlot(mes.epi,group.by = 'SendingType',cols = color_pals$type_colors,raster=F,shuffle = T)
p6 <- DimPlot(mes.epi,group.by = 'ReceivingType',cols = color_pals$type_colors,raster=F,shuffle = T)

plot_grid(p1,p2,p3,p4,p5,p6)

pdf(file='MesenchymalToEpitheliumCommunication_Subset.UMAPS.pdf',width=20,height=10)
print(plot_grid(p1,p2,p3,p4,p5,p6))
dev.off()
## quartz_off_screen 
##                 2

We can also cluster these data, which is rapdily becoming my favorite thing to do:

mes.epi <- FindNeighbors(mes.epi,dims = 1:20)
mes.epi <- FindClusters(mes.epi,res=0.1)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11149
## Number of edges: 369367
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9313
## Number of communities: 4
## Elapsed time: 1 seconds
p1 <- DimPlot(mes.epi,group.by = 'seurat_clusters')

# Take a quick look at how metadata slots compare
table(mes.epi$seurat_clusters,mes.epi$Dataset2.Sending)
##    
##     Quad_E Tri_L
##   0   1380  5340
##   1   2308    42
##   2    876  1063
##   3     81    59
table(mes.epi$seurat_clusters,mes.epi$orig.ident) # We probably want to take a look at this in a plot
##    
##     BEF12 BEF14 BEF15 BEFM1 BEFM2
##   0  1030  2230  2080    63  1317
##   1     8    13    21  2298    10
##   2   344   299   420   121   755
##   3     4    36    19    29    52
table(mes.epi$SendingType,mes.epi$orig.ident) # Note: No pericyte signaling captured in BEFM1...why? Good thing / Bad thing? Is this insightful/useful information?
##                        
##                         BEF12 BEF14 BEF15 BEFM1 BEFM2
##   Aberrant_Fibroblast     350   806   620   165   788
##   Alveolar_Fibroblast     350   806   781  1502   788
##   Cycling_Mes             350   651   730   769   214
##   Pericytes               116    91    50     0   176
##   Remodeling_Fibroblast   220   224   359    75   168
# This measures how different archetypes are distributed over each tissue sample
to.plot <- data.frame(table(mes.epi$seurat_clusters,mes.epi$orig.ident))
p2 <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  #geom_bar(stat='identity',position = 'dodge')+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5)+
  ggtitle('NICHES Archetype Distribution Over Sample')

# This measures how different sending cell populations are distributed over each tissue sample
p3 <- DimPlot(mes.epi,group.by = 'SendingType',cols = color_pals$type_colors,raster=F,shuffle = T)
to.plot <- data.frame(table(mes.epi$SendingType,mes.epi$orig.ident))
# Order the factors for plotting
temp <- to.plot %>% group_by(Var1) %>% summarise(Frequency = sum(Freq))
temp <- arrange(temp,desc(Frequency))
to.plot$Var1 <- factor(to.plot$Var1,levels = temp$Var1)
p4 <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5,)+
  scale_fill_manual(values=color_pals$type_colors)+
  ggtitle('SendingType Distribution Over Sample')

# This measures how different receiving cell populations are distributed over each tissue sample
p5 <- DimPlot(mes.epi,group.by = 'ReceivingType',cols = color_pals$type_colors,raster=F,shuffle = T)
to.plot <- data.frame(table(mes.epi$ReceivingType,mes.epi$orig.ident))
# Order the factors for plotting
temp <- to.plot %>% group_by(Var1) %>% summarise(Frequency = sum(Freq))
temp <- arrange(temp,desc(Frequency))
to.plot$Var1 <- factor(to.plot$Var1,levels = temp$Var1)
p6 <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5,)+
  scale_fill_manual(values=color_pals$type_colors)+
  ggtitle('ReceivingType Distribution Over Sample')

plot_grid(p1,p3,p5,p2,p4,p6)

pdf(file='MesenchymalToEpitheliumCommunication_Subset_Clustered_Distributions.pdf',width=20,height=10)
plot_grid(p1,p3,p5,p2,p4,p6)
dev.off()
## quartz_off_screen 
##                 2

I think this is a cool new plotting technique. What do you think?

This is very cool! Looks like there is a clear effect on Mesenchymal-Epithelial communication due to the addition of Macrophages. Let’s dig into this specific signal a bit by running a statistical test comparing these two categories.

Couple things here. 1. There is a clear association between alveolar fibroblasts with Quad culture, that is already known I think from your analysis but it is cool to see here. 2. Seems like aberrant fibroblasts go down with the addition of macrophages, which again is known from your analysis but i also think is cool. 3. There is a small but possibly real increase in BEFM1 of signaling landing on BASCs. What does this mean / how do we interpret? 4. The BEFM1 archetype seems to be deeply correlated with Alveolar_Fibroblast development. 5. I think that signaling archetype 0 is associated with the Aberrant_Fibroblast phenotype. Let’s check:

# This measures how different sending cell populations are distributed over each signaling archetype
to.plot <- data.frame(table(mes.epi$SendingType,mes.epi$seurat_clusters))
# Order the factors for plotting
temp <- to.plot %>% group_by(Var1) %>% summarise(Frequency = sum(Freq))
temp <- arrange(temp,desc(Frequency))
to.plot$Var1 <- factor(to.plot$Var1,levels = temp$Var1)

plot <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5,)+
  scale_fill_manual(values=color_pals$type_colors)+
  ggtitle('SendingType Distribution Over Archetype')+
  xlab('Signaling Archetype')
plot_grid(p4,plot,nrow=1)

Huh. Not what I was expecting. I think instead, what we are seeing is that Archetype 1 preferentially does not include Aberrant_Fibroblasts (Aberrant_Fibroblasts do not generally communicate with epithelium using the mechanisms present in signaling Archetype 1.) We already know that Archetype 1 is associated with BEFM1. Very cool!!

There are a couple of ways that we can parse out specific signals that we might find to be important biologically. We can leverage any meta-data slot(s) that we want to perform these calculations We could directly compare Tri Culture vs. Quad Culture within this data set, as follows:

Compare MesEpi with and without macrophages:

First, we want a new meta data category that annotates Tri culture vs. Quad culture. we do this as follows:

mes.epi$Dataset4 <- NA
Idents(mes.epi) <- mes.epi$Dataset2.Sending 
mes.epi <- RenameIdents(mes.epi,
                        'Tri_E'='Tri',
                        'Tri_L'='Tri',
                        'Quad_E'='Quad',
                        'Quad_L'='Quad')
mes.epi$Dataset4 <- Idents(mes.epi)
table(mes.epi$Dataset4)
## 
##  Tri Quad 
## 6504 4645

Now, we can run a marker test comparing these two categories specifically:

mark.mes.epi <- FindAllMarkers(mes.epi,only.pos = T,min.pct = 0.1,logfc.threshold = 0.01)
mark.mes.epi$ratio <- mark.mes.epi$pct.1/mark.mes.epi$pct.2
mark.mes.epi$power <- mark.mes.epi$ratio*mark.mes.epi$avg_log2FC
kable(mark.mes.epi %>% group_by(cluster) %>% top_n(50,ratio) %>% arrange(desc(ratio)))
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene ratio power
0 0.0406835 0.255 0.000 0 Quad Wnt7a—Lrp6 Inf Inf
0 0.0359223 0.117 0.001 0 Quad Sct—Sctr 117.000000 4.2029065
0 0.3779125 0.558 0.014 0 Quad Adipoq—Adipor1 39.857143 15.0625137
0 0.0848352 0.209 0.006 0 Quad Lamb3—Itga2 34.833333 2.9550926
0 0.2268704 0.326 0.012 0 Quad Adipoq—Adipor2 27.166667 6.1633118
0 0.0766549 0.225 0.010 0 Quad Lamc2—Itga2 22.500000 1.7247360
0 0.0433286 0.194 0.010 0 Quad Lama1—Itga2 19.400000 0.8405755
0 0.0257216 0.167 0.009 0 Quad Tnfsf12—Tnfrsf25 18.555556 0.4772778
0 0.0279528 0.230 0.013 0 Quad Csf2—Il3ra 17.692308 0.4945493
0 0.0598572 0.192 0.012 0 Quad Cdh1—Egfr 16.000000 0.9577156
0 0.2540416 0.374 0.024 0 Quad Lamb3—Itga3 15.583333 3.9588143
0 0.0269311 0.139 0.009 0 Quad Csf2—Csf2ra 15.444444 0.4159356
0 0.0452318 0.126 0.009 0 Tri Rtn4—Tnfrsf19 14.000000 0.6332454
0 0.0724622 0.206 0.015 0 Quad S100a9—Tlr4 13.733333 0.9951475
0 0.0425042 0.143 0.012 0 Quad Cdh1—Igf1r 11.916667 0.5065085
0 0.1358085 0.357 0.030 0 Quad Lama1—Itga3 11.900000 1.6161206
0 0.2117413 0.373 0.033 0 Tri Gdf10—Acvr1b 11.303030 2.3933186
0 0.2256925 0.396 0.036 0 Quad Lamc2—Itga3 11.000000 2.4826174
0 0.0253900 0.198 0.018 0 Quad Fgf1—Egfr 11.000000 0.2792896
0 0.0202926 0.110 0.010 0 Quad Tcn2—Lrp2 11.000000 0.2232181
0 0.1194167 0.247 0.023 0 Quad Lamb3—Itga6 10.739130 1.2824314
0 0.1925922 0.349 0.033 0 Quad Lamb3—Itgb4 10.575758 2.0368086
0 0.0657472 0.116 0.011 0 Quad Col5a3—Sdc3 10.545454 0.6933345
0 0.0538167 0.175 0.017 0 Quad Wnt7b—Lrp5 10.294118 0.5539958
0 0.3899451 0.457 0.045 0 Quad Lamb3—Itgb1 10.155556 3.9601089
0 0.0262222 0.101 0.010 0 Quad Mdk—Lrp2 10.100000 0.2648437
0 0.1282287 0.371 0.038 0 Quad Lama1—Itgb4 9.763158 1.2519167
0 0.0466438 0.107 0.011 0 Quad Apoe—Lrp2 9.727273 0.4537166
0 0.1443302 0.435 0.046 0 Quad Lama1—Sdc4 9.456522 1.3648620
0 0.0668396 0.217 0.023 0 Quad Cdh1—Lrp5 9.434783 0.6306173
0 0.0565876 0.113 0.012 0 Quad Mdk—Sdc3 9.416667 0.5328662
0 0.0128408 0.129 0.014 0 Quad Fgf1—Fgfr3 9.214286 0.1183191
0 0.0857056 0.284 0.031 0 Quad Angptl4—Tie1 9.161290 0.7851740
0 0.0901725 0.279 0.031 0 Quad Lama1—Itga6 9.000000 0.8115523
0 0.1540670 0.301 0.034 0 Quad Csf2—Itgb1 8.852941 1.3639458
0 0.2538694 0.494 0.058 0 Quad Lama1—Itgb1 8.517241 2.1622666
0 0.0529281 0.110 0.013 0 Quad Wnt5a—Mcam 8.461538 0.4478533
0 0.1192592 0.253 0.030 0 Quad Cdh1—Ptprf 8.433333 1.0057523
0 0.0723921 0.124 0.015 0 Quad Col5a1—Sdc3 8.266667 0.5984410
0 0.1745612 0.179 0.022 0 Tri Slit1—Sdc1 8.136364 1.4202936
0 0.0319312 0.160 0.020 0 Quad Psen1—Notch4 8.000000 0.2554493
0 0.1649228 0.367 0.046 0 Quad Lamc2—Itgb4 7.978261 1.3157968
0 0.1016286 0.261 0.033 0 Quad Lamc2—Itga6 7.909091 0.8037901
0 0.3426242 0.479 0.061 0 Quad Lamc2—Itgb1 7.852459 2.6904429
0 0.0823475 0.101 0.013 0 Quad Anxa1—Dysf 7.769231 0.6397771
0 0.0524192 0.246 0.032 0 Quad Adam15—Itga5 7.687500 0.4029727
0 0.0837317 0.160 0.021 0 Quad Tgfb1—Cxcr4 7.619048 0.6379560
0 0.0227634 0.107 0.015 0 Quad Adm—Calcrl 7.133333 0.1623789
0 0.0461006 0.260 0.037 0 Quad Fgf1—Cd44 7.027027 0.3239505
0 0.2001768 0.219 0.032 0 Quad Mmp9—Ephb2 6.843750 1.3699600
0 0.0230529 0.116 0.017 0 Quad Sema3f—Nrp1 6.823529 0.1573019
0 0.0513693 0.199 0.030 0 Quad Cdh1—Erbb3 6.633333 0.3407498
0 0.0170333 0.131 0.021 0 Quad Lama1—Gpc1 6.238095 0.1062552
0 0.0322468 0.110 0.026 0 Tri Thbs1—Itga2b 4.230769 0.1364290
0 0.1465410 0.139 0.033 0 Tri Fn1—Itga2b 4.212121 0.6172484
0 0.1324566 0.139 0.033 0 Tri Calr—Itga2b 4.212121 0.5579234
0 0.0750097 0.168 0.043 0 Tri Hdc—Hrh4 3.906977 0.2930610
0 0.0461880 0.206 0.054 0 Tri Slit3—Robo2 3.814815 0.1761988
0 0.0516654 0.247 0.066 0 Tri Slit2—Robo2 3.742424 0.1933539
0 0.3357056 0.228 0.067 0 Tri Fn1—Nt5e 3.402985 1.1424013
0 0.2020212 0.345 0.104 0 Tri Col14a1—Cd44 3.317308 0.6701666
0 0.0358461 0.111 0.034 0 Tri Pigf—Flt1 3.264706 0.1170270
0 0.0320742 0.127 0.039 0 Tri Ngf—Ngfr 3.256410 0.1044467
0 0.1303773 0.189 0.059 0 Tri B2m—Cd247 3.203390 0.4176493
0 0.0578606 0.139 0.044 0 Tri Gnai2—Ptpru 3.159091 0.1827870
0 0.0516398 0.147 0.051 0 Tri Tnc—Nt5e 2.882353 0.1488441
0 0.1464161 0.335 0.127 0 Tri Bdnf—Ngfr 2.637795 0.3862157
0 0.2108339 0.129 0.050 0 Tri Fn1—Itga8 2.580000 0.5439515
0 0.1786970 0.361 0.140 0 Tri App—Ngfr 2.578571 0.4607830
0 0.1963006 0.360 0.140 0 Tri Rtn4—Ngfr 2.571429 0.5047729
0 0.6967932 0.781 0.318 0 Tri Thbs2—Itgb1 2.455975 1.7113067
0 0.5210727 0.732 0.301 0 Tri Thbs2—Cd47 2.431894 1.2671935
0 0.0356187 0.245 0.104 0 Tri Mst1—Mst1r 2.355769 0.0839094
0 0.0710921 0.254 0.110 0 Tri Efemp1—Egfr 2.309091 0.1641582
0 0.0902954 0.129 0.056 0 Tri A2m—Lrp1 2.303571 0.2080018
0 0.4032791 0.414 0.183 0 Tri Thbs2—Itga6 2.262295 0.9123364
0 0.1230577 0.166 0.079 0 Tri Dlk1—Notch1 2.101266 0.2585769
0 0.0882803 0.182 0.088 0 Tri Dlk1—Notch2 2.068182 0.1825798
0 0.0172538 0.153 0.075 0 Tri Bmp8b—Bmpr2 2.040000 0.0351978
0 0.1340253 0.370 0.182 0 Tri Icam4—Itgb1 2.032967 0.2724691
0 0.2346354 0.370 0.189 0 Tri Fn1—Itgb7 1.957672 0.4593391
0 0.0166589 0.174 0.096 0 Tri Bmp8b—Acvr2a 1.812500 0.0301943
0 0.0221861 0.175 0.099 0 Tri Bmp8b—Bmpr1a 1.767677 0.0392179
0 0.2945381 0.512 0.290 0 Tri Sfrp1—Fzd6 1.765517 0.5200120
0 0.0117863 0.152 0.087 0 Tri Bmp8b—Acvr1 1.747126 0.0205921
0 0.1396410 0.255 0.147 0 Tri Efemp2—Lingo1 1.734694 0.2422344
0 0.0260303 0.109 0.063 0 Tri Apoe—Lrp8 1.730159 0.0450365
0 0.7234671 0.818 0.473 0 Tri Ptn—Sdc1 1.729387 1.2511545
0 0.0966098 0.254 0.147 0 Tri Rtn4—Lingo1 1.727891 0.1669312
0 0.6257268 0.762 0.442 0 Tri Ptn—Plxnb2 1.723982 1.0787417
0 0.0619612 0.517 0.302 0 Tri Rgma—Bmpr2 1.711920 0.1060726
0 0.1020057 0.241 0.141 0 Tri Vcam1—Itgb1 1.709220 0.1743501
0 0.1064315 0.545 0.327 0 Tri Bmp5—Bmpr2 1.666667 0.1773859
0 0.0132838 0.125 0.075 0 Tri Dkk2—Kremen2 1.666667 0.0221396
0 0.0640236 0.124 0.078 0 Tri Fn1—Itgb8 1.589744 0.1017811
0 0.0234753 0.124 0.078 0 Tri Col4a1—Itgb8 1.589744 0.0373197
0 0.1158296 0.335 0.213 0 Tri Gnai2—S1pr5 1.572770 0.1821733
0 0.0591788 0.260 0.166 0 Tri Gnai2—Adcy9 1.566265 0.0926896
0 0.2718472 0.655 0.424 0 Tri Fbln1—Itgb1 1.544811 0.4199527
0 0.1321240 0.598 0.389 0 Tri Bmp5—Bmpr1a 1.537275 0.2031109

Not bad! Some very interesting signals in there. Let’s remove signaling mediated by integrins, which are less interesting for this project.

MOI <- mark.mes.epi[-grep('Itg',mark.mes.epi$gene),]
MOI <- MOI %>% group_by(cluster) %>% top_n(50,ratio) %>% arrange(desc(ratio))
kable(MOI)
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene ratio power
0 0.0406835 0.255 0.000 0 Quad Wnt7a—Lrp6 Inf Inf
0 0.0359223 0.117 0.001 0 Quad Sct—Sctr 117.000000 4.2029065
0 0.3779125 0.558 0.014 0 Quad Adipoq—Adipor1 39.857143 15.0625137
0 0.2268704 0.326 0.012 0 Quad Adipoq—Adipor2 27.166667 6.1633118
0 0.0257216 0.167 0.009 0 Quad Tnfsf12—Tnfrsf25 18.555556 0.4772778
0 0.0279528 0.230 0.013 0 Quad Csf2—Il3ra 17.692308 0.4945493
0 0.0598572 0.192 0.012 0 Quad Cdh1—Egfr 16.000000 0.9577156
0 0.0269311 0.139 0.009 0 Quad Csf2—Csf2ra 15.444444 0.4159356
0 0.0452318 0.126 0.009 0 Tri Rtn4—Tnfrsf19 14.000000 0.6332454
0 0.0724622 0.206 0.015 0 Quad S100a9—Tlr4 13.733333 0.9951475
0 0.0425042 0.143 0.012 0 Quad Cdh1—Igf1r 11.916667 0.5065085
0 0.2117413 0.373 0.033 0 Tri Gdf10—Acvr1b 11.303030 2.3933186
0 0.0253900 0.198 0.018 0 Quad Fgf1—Egfr 11.000000 0.2792896
0 0.0202926 0.110 0.010 0 Quad Tcn2—Lrp2 11.000000 0.2232181
0 0.0657472 0.116 0.011 0 Quad Col5a3—Sdc3 10.545454 0.6933345
0 0.0538167 0.175 0.017 0 Quad Wnt7b—Lrp5 10.294118 0.5539958
0 0.0262222 0.101 0.010 0 Quad Mdk—Lrp2 10.100000 0.2648437
0 0.0466438 0.107 0.011 0 Quad Apoe—Lrp2 9.727273 0.4537166
0 0.1443302 0.435 0.046 0 Quad Lama1—Sdc4 9.456522 1.3648620
0 0.0668396 0.217 0.023 0 Quad Cdh1—Lrp5 9.434783 0.6306173
0 0.0565876 0.113 0.012 0 Quad Mdk—Sdc3 9.416667 0.5328662
0 0.0128408 0.129 0.014 0 Quad Fgf1—Fgfr3 9.214286 0.1183191
0 0.0857056 0.284 0.031 0 Quad Angptl4—Tie1 9.161290 0.7851740
0 0.0529281 0.110 0.013 0 Quad Wnt5a—Mcam 8.461538 0.4478533
0 0.1192592 0.253 0.030 0 Quad Cdh1—Ptprf 8.433333 1.0057523
0 0.0723921 0.124 0.015 0 Quad Col5a1—Sdc3 8.266667 0.5984410
0 0.1745612 0.179 0.022 0 Tri Slit1—Sdc1 8.136364 1.4202936
0 0.0319312 0.160 0.020 0 Quad Psen1—Notch4 8.000000 0.2554493
0 0.0823475 0.101 0.013 0 Quad Anxa1—Dysf 7.769231 0.6397771
0 0.0837317 0.160 0.021 0 Quad Tgfb1—Cxcr4 7.619048 0.6379560
0 0.0227634 0.107 0.015 0 Quad Adm—Calcrl 7.133333 0.1623789
0 0.0461006 0.260 0.037 0 Quad Fgf1—Cd44 7.027027 0.3239505
0 0.2001768 0.219 0.032 0 Quad Mmp9—Ephb2 6.843750 1.3699600
0 0.0230529 0.116 0.017 0 Quad Sema3f—Nrp1 6.823529 0.1573019
0 0.0513693 0.199 0.030 0 Quad Cdh1—Erbb3 6.633333 0.3407498
0 0.0170333 0.131 0.021 0 Quad Lama1—Gpc1 6.238095 0.1062552
0 0.0349009 0.111 0.018 0 Quad Serpine1—Lrp2 6.166667 0.2152221
0 0.0141602 0.117 0.019 0 Quad Bmp3—Acvr1 6.157895 0.0871973
0 0.0295976 0.134 0.022 0 Quad Bmp3—Bmpr1a 6.090909 0.1802765
0 0.0846650 0.224 0.038 0 Quad S100a8—Tlr4 5.894737 0.4990782
0 0.0242458 0.127 0.022 0 Quad Bmp3—Acvr2a 5.772727 0.1399646
0 0.0221524 0.114 0.020 0 Quad Lrpap1—Lrp2 5.700000 0.1262684
0 0.0537176 0.296 0.053 0 Quad Angpt1—Tie1 5.584906 0.3000077
0 0.0880642 0.111 0.020 0 Quad Ptn—Sdc3 5.550000 0.4887564
0 0.0157701 0.100 0.019 0 Quad Bmp3—Bmpr2 5.263158 0.0830005
0 0.0202024 0.146 0.028 0 Quad Angptl2—Tie1 5.214286 0.1053409
0 0.0980915 0.191 0.040 0 Quad Mmp9—Lrp1 4.775000 0.4683869
0 0.1495824 0.343 0.072 0 Quad Tgfb1—Eng 4.763889 0.7125938
0 0.0145269 0.128 0.027 0 Quad Fgf1—Fgfr2 4.740741 0.0688684
0 0.2366882 0.315 0.067 0 Quad Mmp9—Cd44 4.701492 1.1127878
0 0.0628092 0.233 0.053 0 Quad Lin7c—Abca1 4.396226 0.2761235
0 0.2203756 0.366 0.085 0 Quad Hbegf—Egfr 4.305882 0.9489113
0 0.0714967 0.262 0.061 0 Quad Ntn1—Unc5b 4.295082 0.3070844
0 0.0750097 0.168 0.043 0 Tri Hdc—Hrh4 3.906977 0.2930610
0 0.0461880 0.206 0.054 0 Tri Slit3—Robo2 3.814815 0.1761988
0 0.0516654 0.247 0.066 0 Tri Slit2—Robo2 3.742424 0.1933539
0 0.3357056 0.228 0.067 0 Tri Fn1—Nt5e 3.402985 1.1424013
0 0.2020212 0.345 0.104 0 Tri Col14a1—Cd44 3.317308 0.6701666
0 0.0358461 0.111 0.034 0 Tri Pigf—Flt1 3.264706 0.1170270
0 0.0320742 0.127 0.039 0 Tri Ngf—Ngfr 3.256410 0.1044467
0 0.1303773 0.189 0.059 0 Tri B2m—Cd247 3.203390 0.4176493
0 0.0578606 0.139 0.044 0 Tri Gnai2—Ptpru 3.159091 0.1827870
0 0.0516398 0.147 0.051 0 Tri Tnc—Nt5e 2.882353 0.1488441
0 0.1464161 0.335 0.127 0 Tri Bdnf—Ngfr 2.637795 0.3862157
0 0.1786970 0.361 0.140 0 Tri App—Ngfr 2.578571 0.4607830
0 0.1963006 0.360 0.140 0 Tri Rtn4—Ngfr 2.571429 0.5047729
0 0.5210727 0.732 0.301 0 Tri Thbs2—Cd47 2.431894 1.2671935
0 0.0356187 0.245 0.104 0 Tri Mst1—Mst1r 2.355769 0.0839094
0 0.0710921 0.254 0.110 0 Tri Efemp1—Egfr 2.309091 0.1641582
0 0.0902954 0.129 0.056 0 Tri A2m—Lrp1 2.303571 0.2080018
0 0.1230577 0.166 0.079 0 Tri Dlk1—Notch1 2.101266 0.2585769
0 0.0882803 0.182 0.088 0 Tri Dlk1—Notch2 2.068182 0.1825798
0 0.0172538 0.153 0.075 0 Tri Bmp8b—Bmpr2 2.040000 0.0351978
0 0.0166589 0.174 0.096 0 Tri Bmp8b—Acvr2a 1.812500 0.0301943
0 0.0221861 0.175 0.099 0 Tri Bmp8b—Bmpr1a 1.767677 0.0392179
0 0.2945381 0.512 0.290 0 Tri Sfrp1—Fzd6 1.765517 0.5200120
0 0.0117863 0.152 0.087 0 Tri Bmp8b—Acvr1 1.747126 0.0205921
0 0.1396410 0.255 0.147 0 Tri Efemp2—Lingo1 1.734694 0.2422344
0 0.0260303 0.109 0.063 0 Tri Apoe—Lrp8 1.730159 0.0450365
0 0.7234671 0.818 0.473 0 Tri Ptn—Sdc1 1.729387 1.2511545
0 0.0966098 0.254 0.147 0 Tri Rtn4—Lingo1 1.727891 0.1669312
0 0.6257268 0.762 0.442 0 Tri Ptn—Plxnb2 1.723982 1.0787417
0 0.0619612 0.517 0.302 0 Tri Rgma—Bmpr2 1.711920 0.1060726
0 0.1064315 0.545 0.327 0 Tri Bmp5—Bmpr2 1.666667 0.1773859
0 0.0132838 0.125 0.075 0 Tri Dkk2—Kremen2 1.666667 0.0221396
0 0.1158296 0.335 0.213 0 Tri Gnai2—S1pr5 1.572770 0.1821733
0 0.0591788 0.260 0.166 0 Tri Gnai2—Adcy9 1.566265 0.0926896
0 0.1321240 0.598 0.389 0 Tri Bmp5—Bmpr1a 1.537275 0.2031109
0 0.3591713 0.585 0.384 0 Tri Gpc3—Cd81 1.523438 0.5471751
0 0.1389964 0.364 0.241 0 Tri Fgf10—Fgfr2 1.510373 0.2099365
0 0.0275825 0.221 0.147 0 Tri Gnai2—Cxcr2 1.503401 0.0414676
0 0.0779625 0.199 0.134 0 Tri Dcn—Met 1.485075 0.1157801
0 0.0802552 0.601 0.405 0 Tri Bmp5—Acvr2a 1.483951 0.1190947
0 0.0402334 0.293 0.200 0 Tri Bmp5—Acvr2b 1.465000 0.0589419
0 0.1508341 0.361 0.248 0 Tri Apoe—Ldlr 1.455645 0.2195609
0 0.0754248 0.523 0.360 0 Tri Bmp5—Acvr1 1.452778 0.1095755
0 0.0171711 0.193 0.134 0 Tri Cntf—Il6r 1.440298 0.0247316
0 0.0694165 0.218 0.153 0 Tri B2m—Cd3g 1.424837 0.0989072
0 0.0560992 0.190 0.136 0 Tri Cxcl12—Ackr3 1.397059 0.0783739
0 0.1617107 0.433 0.314 0 Tri Cthrc1—Fzd6 1.378981 0.2229959

So there is definitely a huge activation of WNT family signaling here. Wnt7a, Wnt7b, via Lrp5/6. Rspo1-Lgr4 (I know the ligand well, but not the receptor!) Also mild upregulation of Fgf1/2/7 signaling via Fgfr3. Also a big upregulation of Hbegf-Egfr. The tri-culture elevation of Dlk1-Notch1 is interesting, not sure what it means. Nrg1-Erbb3 also upregulated in Quad, and Erbb3 is an powerful epithelial stem-cell / growth activator. The Csf2-Csf2ra upregulated in Quad might be relevant to macrophage biology. Efna4—Epha7 also upregulated, it is a spatial guidance molecule implicated in tissue morphogenesis, as is Sema4b—Dcbld2, which is acting through a receptor I have never seen before. Il11—Il6st also upregulated. Immune recruitment / activation? Not sure what this mechanism does, would need to read.

Let’s make some feature plots of Mes-Epi signals that are upregulated in Quad culture as compared to Tri, like this:

quick.function <- function(mech){
p1 <- FeaturePlot(mes.epi,mech)
p2 <- VlnPlot(mes.epi,mech)
p3 <- VlnPlot(mes.epi,mech,group.by = 'orig.ident')
plot_grid(p1,p2,p3,nrow = 1)
}
mech <- 'Fgf7—Fgfr3'
quick.function(mech)

mech <- 'Wnt7a—Lrp6'
quick.function(mech)

mech <- 'Wnt7b—Lrp5'
quick.function(mech)

mech <- 'Csf2—Csf2ra'
quick.function(mech)

mech <- 'Fgf2—Fgfr3'
quick.function(mech)

mech <- 'Fgf2—Fgfr2'
quick.function(mech)

mech <- 'Hbegf—Egfr'
quick.function(mech)

mech <- 'Il11—Il6st'
quick.function(mech)

mech <- 'Rspo1—Lgr4'
quick.function(mech)

EDITS 11/9/2023, post-subset: I have added a third plot here that looks at the data split by Tissue, so that we can see if any of the effects are sample-specific. So, while these mechanisms are definitely still interesting, I am seeing the same thing you have mentioned, which is that BEFM1 is like “much better” than other quad cultures. Weird. But I guess also makes sense? We already kind of knew this. Not a huge issue, the tissue looks the best too…

Next, let’s look at some Tri culture Mes-Epi markers:

Not sure what Dlk1 is doing or if it is relevant, it is a Notch ligand and it is definitely lowered in BEFM1:

mech <- 'Dlk1—Notch1'
quick.function(mech)

Thbs2 is Themis’ pro-fibrotic, anti-angiogenic matrix molecule that he knocks out, wild that it is here so clearly and is so reduced in BEFM1 and the quad cultures more generally

mech <- 'Thbs2—Cd47'
quick.function(mech)

Sfrp1—Fzd6 is an extracellular binder of WNT ligands that sequesters them outside of cells and rpevents them from binding, it is net anti-WNT. It is DOWNREGULATED in quad culture, suggesting that macrophage addition will activate this pathway and promote epithelial WNT signal transduction

mech <-'Sfrp1—Fzd6'
quick.function(mech)

Cthrc1—Fzd6 is a major matrix molecule implicated in fibrosis (ligand): https://www.nature.com/articles/s41467-020-15647-5 and the receptor is a critical cue for alveolar regeneration by epithelium: https://www.cell.com/cell/pdf/S0092-8674(23)00541-X.pdf

mech <-'Cthrc1—Fzd6'
quick.function(mech)

But I prefer to look at how the signaling archetypes difer from one another. It is more generalizable, and also it leverages the actual data structure to perform the entire analysis, which means that the signals that come out as output mechanisms are much clearer observed in the data itself. It is easier to interpret this way. So, here is a parallel workflow which uses the archetype information as the main data handle:

Differential testing by Archetype

Idents(mes.epi) <- mes.epi$seurat_clusters # This makes the archetype information the main handle
mark.mes.epi <- FindAllMarkers(mes.epi,only.pos = T,min.pct = 0.1,logfc.threshold = 0.01)
mark.mes.epi$ratio <- mark.mes.epi$pct.1/mark.mes.epi$pct.2
mark.mes.epi$power <- mark.mes.epi$ratio*mark.mes.epi$avg_log2FC
kable(mark.mes.epi %>% group_by(cluster) %>% top_n(10,ratio))
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene ratio power
0 0.3664661 0.391 0.084 0 0 App—Ngfr 4.654762 1.7058126
0 0.3373300 0.390 0.084 0 0 Rtn4—Ngfr 4.642857 1.5661749
0 0.2120839 0.361 0.078 0 0 Bdnf—Ngfr 4.628205 0.9815679
0 0.0599403 0.278 0.048 0 0 Mst1—Mst1r 5.791667 0.3471544
0 0.1153768 0.278 0.057 0 0 Gnai2—Cxcr2 4.877193 0.5627149
0 0.0899865 0.221 0.047 0 0 Cxcl1—Cxcr2 4.702128 0.4231281
0 0.0399730 0.215 0.049 0 0 Vasp—Cxcr2 4.387755 0.1753916
0 0.0743570 0.150 0.034 0 0 Dkk2—Kremen2 4.411765 0.3280457
0 0.0379441 0.129 0.032 0 0 Ngf—Ngfr 4.031250 0.1529620
0 0.0262909 0.110 0.025 0 0 Vcam1—Itgb7 4.400000 0.1156799
0 0.6116117 0.941 0.054 0 1 Adipoq—Adipor1 17.425926 10.6579010
0 0.4329424 0.651 0.041 0 1 Lamb3—Itga3 15.878049 6.8742809
0 0.3752132 0.541 0.036 0 1 Adipoq—Adipor2 15.027778 5.6386208
0 0.1516485 0.373 0.015 0 1 Lamb3—Itga2 24.866667 3.7709915
0 0.1404703 0.409 0.017 0 1 Lamc2—Itga2 24.058823 3.3795490
0 0.0736566 0.475 0.008 0 1 Wnt7a—Lrp6 59.375000 4.3733602
0 0.0681361 0.220 0.003 0 1 Sct—Sctr 73.333333 4.9966447
0 0.0470222 0.241 0.016 0 1 Csf2—Csf2ra 15.062500 0.7082716
0 0.0469633 0.356 0.023 0 1 Fgf1—Egfr 15.478261 0.7269106
0 0.0543063 0.120 0.007 0 1 Calcb—Ramp1 17.142857 0.9309652
0 0.1040387 0.294 0.024 0 2 Wnt5a—Ror1 12.250000 1.2744745
0 0.0834049 0.272 0.013 0 2 Cntf—Lifr 20.923077 1.7450864
0 0.0631832 0.238 0.022 0 2 Pon2—Htr2a 10.818182 0.6835278
0 0.2816088 0.157 0.003 0 2 Fn1—Col13a1 52.333333 14.7375259
0 0.1553227 0.205 0.017 0 2 Gnai2—Adra2b 12.058824 1.8730091
0 0.1008910 0.140 0.003 0 2 Nid1—Col13a1 46.666667 4.7082451
0 0.0531654 0.137 0.002 0 2 Nid2—Col13a1 68.500000 3.6418282
0 0.0658717 0.141 0.004 0 2 Icam1—Il2ra 35.250000 2.3219768
0 0.0269528 0.135 0.009 0 2 Ctf1—Lifr 15.000000 0.4042924
0 0.1260287 0.110 0.010 0 2 Col1a1—Itga11 11.000000 1.3863162
0 0.6905106 0.907 0.007 0 3 Vegfa—Itga9 129.571429 89.4704476
0 0.6463646 0.921 0.007 0 3 Spp1—Itga9 131.571429 85.0431074
0 0.3968352 0.693 0.005 0 3 Tnc—Itga9 138.600000 55.0013534
0 0.3476353 0.929 0.007 0 3 Vegfc—Itga9 132.714286 46.1361757
0 0.2938716 0.664 0.004 0 3 Adam12—Itga9 166.000000 48.7826932
0 0.0902893 0.343 0.002 0 3 Adam15—Itga9 171.500000 15.4846157
0 0.0677749 0.186 0.001 0 3 Vcam1—Itga9 186.000000 12.6061392
0 0.0207233 0.136 0.001 0 3 Csf2—Itga9 136.000000 2.8183676
0 0.0125385 0.150 0.001 0 3 Wnt7a—Fzd9 150.000000 1.8807717
0 0.0308627 0.129 0.001 0 3 Lamb3—Cd151 129.000000 3.9812839

When we then look at feature plots, things look very nice. We can view trends in the data that correlate with sample, but we didn’t use those labels to perform the analysis. This keeps the analysis indepedent and non-biased.

FeaturePlot(mes.epi,'Wnt7a—Lrp6',max.cutoff = 0.2,split.by = 'orig.ident')

FeaturePlot(mes.epi,'Csf2—Csf2ra',max.cutoff = 0.5,split.by = 'orig.ident',order = T)

FeaturePlot(mes.epi,'Sct—Sctr',max.cutoff = 0.5,split.by = 'orig.ident',order = T) # What is this molecular interaction? It's wild.

FeaturePlot(mes.epi,'Cxcl1—Cxcr2',max.cutoff = 0.9,split.by = 'orig.ident',order=T) # BEFM1 has reduced inflammation ????

FeaturePlot(mes.epi,'Wnt5a—Ror1',max.cutoff = 0.9,split.by = 'orig.ident',order=T)

Reduced Wnt5a in BEFM ???? I think of Wnt5a as a “pro-differentation” cue…could the BEFM1 condition be promoting stemness in the epithelial cells? Is this possible / supported by your analysis?

Downregulated Wnt5a and upregulated Wnt7a … connected to the increased BASC representation in this archetype?

If we want to, we can create a complex heatmap which shows a great deal of information at once, something like as follows:

Let’s make a heatmap that displays all key findings for Mes-Epi Signaling in this dataset

source('/Users/msbr/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/My Drive/Tuft_Sox9_Pneumonectomy_Project/Homeostatic_Single_Cell/CustomHeatmap.R')

mechs.to.annotate <- mark.mes.epi %>% group_by(cluster) %>% top_n(10,ratio)
CustomHeatmap(object = mes.epi,
              data.type = 'CellToCell',
              primary = 'seurat_clusters' ,
              secondary = 'SendingType' ,
              tertiary = 'ReceivingType' ,
              quarternary = 'orig.ident' ,
              primary.cols = NULL,
              secondary.cols = color_pals$type_colors, # Need to be a named list of colors in the right order
              tertiary.cols = color_pals$type_colors,
              quarternary.cols = color_pals$sample_colors,
              features = unique(mark.mes.epi$gene),
              labels = c('Signaling Archetype','Sending Cell Type','Receiving Cell Type','Tissue'),
              selected.row.anotations = mechs.to.annotate$gene,
              selected.label.size = 10)

As a weird sort of cross-project scientific experiment, what happens if I ask to annotate code-inheritied markers of interest from the pneumonectomy project??

mechs.to.annotate <- c('Il17b—Il17rb','Il25—Il17rb',
                                          'Ereg—Erbb3','Areg—Erbb3',
                                          'Sema4d—Plxnb2','Sema6a—Plxna2',#'Sema6a—Plxna4',
                                          'Wnt3a—Fzd6','Wnt7b—Lrp5','Wnt5a—Fzd5','Wnt5a—Fzd6','Wnt3a—Ryk',
                                          'Rspo2—Lgr5','Rspo3—Sdc4','Rspo1—Lrp6','Rspo1—Lrp6','Rspo1—Lgr5',
                                          'Dll1—Notch2','Dll4—Notch2','Dll3—Notch2',
                                          'Sct—Sctr','Dhh—Ptch1',
                                          'Tgfb1—Tgfbr1','Tgfb1—Tgfbr2','Ccl2—Ackr4','Ccl5—Sdc1',#'Ccl5—Sdc4',
                                          'Mdk—Sdc4','Ccl11—Ackr4','Fbln1—Itgb1',
                                          'Il18—Cd48','Tnf—Ltbr','Igf1—Insr','Tnfsf13—Tnfrsf13b')
CustomHeatmap(object = mes.epi,
              data.type = 'CellToCell',
              primary = 'seurat_clusters' ,
              secondary = 'SendingType' ,
              tertiary = 'ReceivingType' ,
              quarternary = 'orig.ident' ,
              primary.cols = NULL,
              secondary.cols = color_pals$type_colors, # Need to be a named list of colors in the right order
              tertiary.cols = color_pals$type_colors,
              quarternary.cols = color_pals$sample_colors,
              features = unique(mark.mes.epi$gene),
              labels = c('Signaling Archetype','Sending Cell Type','Receiving Cell Type','Tissue'),
              selected.row.anotations = mechs.to.annotate,
              selected.label.size = 10)

Fascinating. Seems like there are some markers that overlap and might be biologcally relevant in both cases. Could we learn from this? Could we adopt a similar apporach to compare engineered biology to native? I need to think on this more. Seems powerful to me though.

My overarching summary, in a nutshell:

Tri-culture has major upregualtion of fibrotic pathways, genes, and matrix molecules. Adding macrophages basically removes this signal.

Adding macrophages in Quad-culture promotes significant upregulation of WNT-signaling and FGF-family signaling. We also see the expression of prominent sptial-guidance molecules crucial to local tissue organziation that we know from normal alveolar homeostasis.

I’m beginning to consider if the BEFM1 condition / Quad-Culture condition / addition of macrophages might promote stemness in these cultures. Would this make sense with the rest of our data?